SciPy

Disclaimer: this course is adapted from the notebooks by

Introduction

SciPy is a scientific library that builds upon NumPy. Among others, SciPy deals with:

Integration (scipy.integrate)
Optimization (scipy.optimize)
Interpolation (scipy.interpolate)
Fourier Transform (scipy.fftpack)
Signal Processing (scipy.signal)
Linear Algebra (scipy.linalg)
Sparse matrices (scipy.sparse)
Statistics (scipy.stats)
Image processing (scipy.ndimage)
IO (input/output) (scipy.io)

%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display

Animation display with `python`

Animation with `matplotlib`

You can use FuncAnimation to animate a sequence of images:

%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")


def init():
    ax.set_xlim(0, 2 * np.pi)
    ax.set_ylim(-1, 1)
    return (ln,)


def update(frame):
    xdata.append(frame)
    ydata.append(np.sin(frame))
    ln.set_data(xdata, ydata)
    return (ln,)

ani = FuncAnimation(
    fig,
    update,
    interval=50,
    frames=np.linspace(0, 2 * np.pi, 100),
    init_func=init,
    blit=True,
)

display(HTML(ani.to_jshtml()))

Another way of displaying video exists, using html5 video:

display(HTML(ani.to_html5_video()))

References:

Matplotlib Animations / JavaScript Widgets by Louis Tiao

Animation with `plotly`

matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:

import plotly.express as px
from plotly.offline import plot

df = px.data.gapminder()
fig = px.scatter(
    df,
    x="gdpPercap",
    y="lifeExp",
    animation_frame="year",
    animation_group="country",
    size="pop",
    color="continent",
    hover_name="country",
    log_x=True,
    size_max=55,
    range_x=[100, 100000],
    range_y=[25, 90],
)
fig.show("notebook")

Linear algebra

scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.

References:

Scipy documentation

Solving linear systems:

Find x such that: A x = b for specified matrix A and vector b.

A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
b = np.array([[1, 2, 3]], dtype=np.float64).T

print(A, b)

x = linalg.solve(A, b)
print(x, x.shape, b.shape)

[[ 1.  0.  3.]
 [ 4.  5. 12.]
 [ 7.  8.  9.]] [[1.]
 [2.]
 [3.]]
[[ 0.8       ]
 [-0.4       ]
 [ 0.06666667]] (3, 1) (3, 1)

Check the result at a given precision (different from ==)

np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)

True

Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!

Eigenvalues/ Eigenvectors

A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:

A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)

np.allclose(A, evecs @ np.diag(evals) @ evecs.T)

[-0.94760646+0.j  2.04440479+0.j  3.47941641+0.j] 
 ------
 [[ 0.50787323 -0.78900444 -0.34574091]
 [-0.85203941 -0.51920701 -0.06673021]
 [ 0.12686067 -0.32847537  0.93595422]]

True

EXERCISE: Eigenvalues/Eigenvectors

Verify numerically that the outputs from linalg.eig are indeed approximately eigenvalues and eigenvectors of matrix A above.

Hint: use Scipy documentation on allclose

Symmetric matrices

If A is symmetric you should use eigvalsh (H for Hermitian) instead: This is more robust and leverages the structures (you know they are real!)

Matrix operations

linalg.trace(A) # trace
linalg.det(A) # determinant
linalg.inv(A) # Inverse, consider NEVER using it though :)

Norms

print(linalg.norm(A, ord="fro"))  # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))

4.145345283920557
4.145345283920557
3.4794164120004525
[0.94760646 2.04440479 3.47941641]
4.144258214937556

EXERCISE: Norms computation

Check numerically what the instruction linalg.norm(A, ord=np.inf) is computing. Double check with the help and a numerical test.

A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))

3.016294053009224

Random generation, distributions, etc.

References:

Good practices with numpy random number generators by Albert Thomas
Numpy documentation on RandomState
Random Widgets, by Joseph Salmon: Visualization of various popular distributions.

seed = 12345
rng = np.random.default_rng(seed)  # can be called without a seed
rng.random()

0.22733602246716966

Optimization

Goal: find functions minima or maxima

References:

Scipy Lectures on mathematical optimization.

from scipy import optimize

Finding (local!) minima

def f(x):
    return 4 * x ** 3 + (x - 2) ** 2 + x ** 4


def mf(x):
    return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)


xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()

Default solver for minimization/maximization: fmin_bfgs (see Wikipedia on BFGS)

x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)


plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()

Optimization terminated successfully.
         Current function value: -3.506641
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8
Optimization terminated successfully.
         Current function value: -6.201654
         Iterations: 5
         Function evaluations: 12
         Gradient evaluations: 6
Optimization terminated successfully.
         Current function value: 2.804988
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8

EXERCISE: Basin of attraction

Draw the points on the curves with two different colors :

orange: for the points leading to find the left local minima
red: for the points leading to the right local minima.

Find the zeros of a function

Find x such that f(x) = 0, with fsolve.

omega_c = 3.0

def f(omega):
    return np.tan(2 * np.pi * omega) - omega_c / omega


x = np.linspace(1e-8, 3.2, 1000)
y = f(x)

# Remove vertical lines when the function flips signs
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)

optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))

my_zeros = (
    np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation

from scipy.optimize import curve_fit


def f(x, a, b, c):
    """f(x) = a exp(-bx) + c."""
    return a * np.exp(-b * x) + c


x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5)  # true signal
yn = y + 0.2 * np.random.randn(len(x))  # noisy added

plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()

(a, b, c), _ = curve_fit(f, x, yn)
print(a,"\n", b,"\n", c)

2.6616069003229623 
 1.4624690685428108 
 0.5538807148955575

Displaying

plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()

For polynomial fitting, one can directly use numpy functionsL

x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()

/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:

Polyfit may be poorly conditioned

Interpolation

from scipy.interpolate import interp1d, CubicSpline


def f(x):
    return np.sin(x)


n = np.arange(0, 10)
x = np.linspace(0, 9, 100)

y_meas = f(n) + 0.1 * np.random.randn(len(n))  # add noise
y_real = f(x)

linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)

cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)


plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()

Images

from scipy import ndimage, datasets

img = datasets.face()
type(img), img.dtype, img.ndim, img.shape


print(2 ** 8)  # uint8-> code sur 256 niveau.

n_1, n_2, n_3 = img.shape
np.unique(img)


# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()

RGB decomposition

fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))

plt.tight_layout()

print(img.flags)  # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags)  # can edit now


img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show()

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : False
  ALIGNED : True
  WRITEBACKIFCOPY : False

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

Convert a color image in grayscale

plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()

Blurring (🇫🇷: floutage)

img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()

Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)

#| eval: false
%matplotlib widget
import ipywidgets as widgets

# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()

@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
    """Remove old lines from plots and plot new ones."""
    # [l.remove() for l in ax.lines]
    img_flou = ndimage.gaussian_filter(img, sigma)
    ax.imshow(img_flou, cmap=plt.cm.gray)

Changing colors in an image

import pooch
import os

url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)

img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)


fig, ax = plt.subplots(3, 2)
fig.set_size_inches(7, 4.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)

plt.tight_layout()

EXERCISE: Make the Brazilian italianer

Create a version of the Brazilian flag as follows:

References:

Introduction to geopandas
Course on images/slides by Joseph Salmon
Official SciPy web page
SciPy User Guide
Scipy lectures
The SciPy source code

Introduction

Animation display with python

Animation with matplotlib

Animation with plotly

Linear algebra

Solving linear systems:

Eigenvalues/ Eigenvectors

Symmetric matrices

Matrix operations

Norms

Random generation, distributions, etc.

Optimization

Finding (local!) minima

Find the zeros of a function

Parameters estimation

Interpolation

Images

RGB decomposition

Convert a color image in grayscale

Blurring (🇫🇷: floutage)

Changing colors in an image

Animation display with `python`

Animation with `matplotlib`

Animation with `plotly`